﻿using MathNet.Numerics;
using SeeSharpTools.JY.DSP.Fundamental;
using SeeSharpTools.JY.DSP.Utility.Filter1D;
using Sunny.UI;
using System;
using System.Linq;

namespace SeeSharpTools.Demo
{
    public partial class FormOnline : UIPage
    {
        /// <summary>
        /// 初始信号
        /// </summary>
        private double[] x;

        /// <summary>
        /// 时间数组
        /// </summary>
        private double[] t;

        /// <summary>
        /// 初始相位，单位为度
        /// </summary>
        private double phase;

        /// <summary>
        /// 滤波器分子系数
        /// </summary>
        private double[] b;

        /// <summary>
        /// 滤波器分母系数
        /// </summary>
        private double[] a;

        private double[] y1;

        private double[] y2;

        private bool isRunning;

        /// <summary>
        /// 滤波器内部状态
        /// </summary>
        private double[] zi;

        private double[] zf;

        public FormOnline()
        {
            InitializeComponent();

            isRunning = false;

            phase = 23; // 初始相位，角度单位
            double sampleRate = 1000;
            t = Generate.LinearRange(0, 1 / sampleRate, 0.065 - 1 / sampleRate);

            (b, a) = IIRDesign.Butter(3, new double[] { 100 }, IIRBandType.Lowpass, sampleRate);
            //zi = new double[b.Length - 1];

            timer1.Interval = 500;
        }

        private void timer1_Tick(object sender, EventArgs e)
        {
            x = Sine(t, 50, 1, ref phase, false);
            double[] noise = new double[x.Length];
            Generation.UniformWhiteNoise(ref noise, 0.1);
            x = x.Zip(noise, (a, b) => a + b).ToArray();

            // 每次循环都用滤波器默认的初始状态，会在边界产生瞬态效应
            y1 = IIRFiltering.IIRFilter(b, a, x);

            // 在数据量较大时可以对其进行分段连续滤波，消除瞬态效应
            // 使用滤波器延迟的初始条件和最终条件对各部分中的数据进行滤波
            // 输出对本次循环进行滤波的最终条件，以便在末尾存储滤波器的内部状态
            // 这次滤波的最终条件作为下一个循环滤波的初始条件
            (y2, zf) = IIRFiltering.IIRFilter(b, a, x, zi);
            zi = zf.ToArray();

            double[,] combined = new double[3, x.Length];
            for (int i = 0; i < x.Length; i++)
            {
                combined[0, i] = x[i];
                combined[1, i] = y1[i];
                combined[2, i] = y2[i];
            }

            stripChartX1.Plot(combined);
            stripChartX1.Series[0].Name = "原始噪声信号";
            stripChartX1.Series[1].Name = "单态模式滤波";
            stripChartX1.Series[2].Name = "连续模式滤波";
        }

        private void uiSymbolButton_Click(object sender, EventArgs e)
        {
            if (isRunning)
            {
                timer1.Stop();
                uiSymbolButton.Symbol = 61515;
                isRunning = false;
            }
            else
            {
                timer1.Start();
                uiSymbolButton.Symbol = 61516;
                isRunning = true;
            }
        }

        /// <summary>
        /// 连续相位正弦波生成
        /// </summary>
        /// <param name="t">采样时间序列</param>
        /// <param name="frequency">正弦波的频率</param>
        /// <param name="amplitude">正弦波的幅度</param>
        /// <param name="phase">正弦波的初始相位，以度为单位</param>
        /// <param name="resetPhase">确定正弦波的初始相位，若为false则输出相位可以作为下一次执行的相位输入，使正弦波相位连续</param>
        /// <returns>输出的正弦波和正弦波下一个采样的相位（单位为度）</returns>
        public static double[] Sine(double[] t, double frequency, double amplitude, ref double phase, bool resetPhase = true)
        {
            double[] x = Sine(t, frequency, amplitude, phase);
            double sampleRate = 1 / (t[1] - t[0]);

            if (!resetPhase)
            {
                // 保留本次输出的相位信息，作为正弦波下一个采样点的相位输入
                double oneperiod = sampleRate / frequency;
                double lastphase = Remainder(x.Length, oneperiod) / oneperiod * 360;
                phase += lastphase;
                phase %= 360;
            }

            return x;
        }

        /// <summary>
        /// 正弦波生成
        /// </summary>
        /// <param name="t">采样时间序列</param>
        /// <param name="frequency">正弦波的频率</param>
        /// <param name="amplitude">正弦波的幅度</param>
        /// <param name="phase">正弦波的初始相位，以度为单位</param>
        /// <returns>输出的正弦波</returns>
        public static double[] Sine(double[] t, double frequency, double amplitude = 1, double phase = 0)
        {
            if (amplitude <= 0) throw new Exception("幅度必须大于0");
            double sampleRate = 1 / (t[1] - t[0]);
            if (frequency > sampleRate / 2) throw new Exception("信号频率不能超过采样率的一半");

            double[] x = new double[t.Length];
            double rad = Trig.DegreeToRadian(phase);
            for (int i = 0; i < x.Length; i++)
            {
                x[i] = amplitude * Math.Sin(2 * Math.PI * frequency * t[i] + rad);
            }
            return x;
        }

        /// <summary>
        /// 计算除后的余数（求余运算），符号和被除数相同
        /// </summary>
        /// <param name="x">被除数</param>
        /// <param name="y">除数</param>
        /// <returns></returns>
        public static double Remainder(double x, double y)
        {
            if (double.IsNaN(x) || double.IsNaN(y)) return double.NaN;
            if (y == 0) return double.NaN;
            if (Math.Abs(x - y) < 1e-10) return 0; // 被除数等于除数

            return x % y;
        }
    }
}